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Abstract.  This  paper  discusses  analytical  and  numerical  issues  related  to  elliptic  equa¬ 
tions  with  random  coefficients  which  are  generally  nonlinear  functions  of  white  noise.  Sin¬ 
gularity  issues  are  avoided  by  using  the  Ito-Skorohod  calculus  to  interpret  the  interactions 
between  the  coefficients  and  the  solution.  The  solution  is  constructed  by  means  of  the 
Wiener  Chaos  (Cameron-Martin)  expansions.  The  existence  and  uniqueness  of  the  solu¬ 
tions  are  established  under  rather  weak  assumptions,  the  main  of  which  requires  only  that 
the  expectation  of  the  highest  order  (differential)  operator  is  a  non-degenerate  elliptic  op¬ 
erator.  The  deterministic  coefficients  of  the  Wiener  Chaos  expansion  of  the  solution  solve  a 
lower-triangular  system  of  linear  elliptic  equations  (the  propagator).  This  structure  of  the 
propagator  insures  linear  complexity  of  the  related  numerical  algorithms.  Using  the  lower 
triangular  structure  and  linearity  of  the  propagator,  the  rate  of  convergence  is  derived  for 
a  spectral/hp  finite  element  approximation.  The  results  of  related  numerical  experiments 
are  presented. 


1.  Introduction 


The  objective  of  this  paper  is  to  study,  both  analytically  and  numerically,  linear  elliptic 
equations  with  random  coefficients.  The  perturbations  introduced  by  the  random  coeffi¬ 
cients  can  represent  actual  physical  input  or  uncertainty  about  the  model.  In  this  paper, 
we  limit  our  considerations  to  coefficients  that  are  (possibly  nonlinear)  functions  of  Gauss¬ 
ian  random  fields.  A  popular  example  of  an  equation  with  random  coefficients  is  Darcy 
equation  for  pressure  with  lognormal  permeability  field  (see  [8]  and  the  references  therein). 
It  can  be  specified  as  the  following  Dirichlet  problem: 

—V  (a  (x,uj)  Vu(x))  =  f(x),  x  G  O  C  Md,  u\qq  =  g  (x)  ,  (1.1) 
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where 


a  (x,  uj) 


a(x)  exp 


-E  |e(*)| 


(1.2) 


uj  symbolizes  “chance”,  e(x)  =  Ylk> iak{x)£k  is  the  noise  term,  Ok{x)  are  deterministic 
functions,  and  £  :=  {£&  (w)}fc>1  is  a  set  of  uncorrelated  Gaussian  random  variables  with  zero 
mean  and  unit  variance.  The  force  /  and  the  boundary  value  g  as  well  as  the  boundary  itself 
are  assumed  to  be  non-random.  Clearly,  e(x)  is  a  Gaussian  random  field  with  zero  mean  and 
covariance  Q(x,y )  =  Y2kak(x)ak(y)  and  a(x)  is  the  mean  of  a(x,uj):  a(x)  =  E[a(aj,o;)]  . 


Another  typical  example  is  (1.1)  with 


a(x,uj)  :=  a(x)  +  e(x,u>)  (1.3) 

Models  (1.1)— (1.3)  have  been  actively  investigated  recently  in  the  context  of  uncertainty 
quantification  problems  (see  e.g  [2]).  Note  that  equation  (1.1)  with  the  function  a(x,uj) 
from  (1.3)  is  ill-posed  because  a(x,uj)  can  take  negative  values. 

Linear  equation  with  multiplicative  noise  are  usually  referred  to  as  bi-linear.  These  equa¬ 
tions  are  much  more  challenging  than  linear  equations  with  additive  random  perturbations. 
In  fact,  equation  (1.1)  is  often  too  singular  to  be  “physical”.  By  this  reason,  certain 
regularization  is  in  order  for  each  of  these  two  models. 

One  effective  and  popular  regularization  procedure  for  these  models  is  to  replace  equation 
by  its  Wick  product  version 

—V  (a  (x,  uj)  o  Vu(x))  =  f(x),  x  e  O  C  Md,  u\qq  =  g  (x) ;  (1.4) 

see  [8]  for  a  discussion  of  the  Wick  product  as  a  modeling  tool.  Admittedly,  “Wick  product” 
is  anything  but  product,  it  is  rather  a  stochastic  convolution.  In  fact,  Wick  product  is  a 
version  of  Malliavin  divergence  operator.  For  example,  a  (x,  uj)  o  X7u(x)  is  the  action  of  the 
Malliavin  divergence  operator  with  respect  to  a(x,uj)  on  the  random  field  \7u(x). 

The  idea  of  this  regularization  procedure  was  championed  by  K.  Ito  in  his  seminal  work  [9]. 
Specifically,  Ito  has  replaced  the  “product”  model 

u(t)  =  a  ( u(t ))  +  b  (u  (t))  ■  W  (■ t ) 


by  a  “stochastic  convolution”  model 


u(t)  =  a  (u(t))  +  b  (u  (t))  o  W  (t) 
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which  turned  out  to  be  equivalent  to  the  stochastic  Ito  differential: 

b  (u  (t))  o  Wdt  =  b  (u  (t))  dW  (t)  . 


In  this  paper  we  extend  Ito’s  approach  to  elliptic  stochastic  PDEs.  Since  absence  of  the 
time  variable  prevents  elliptic  SPDEs  from  being  “causal”,  the  underlying  Ito  integrals  are 
replaced  by  the  Skorohod  integrals  (see  [15]).  The  subsequent  analysis  of  the  equation  is 
based  on  the  Wiener  chaos  expansion  (WCE)  with  respect  to  the  Cameron- Martin  basis. 
The  Cameron-Martin  basis  consists  of  random  variables 


where  a  =  (au,  ctQ;  ■  •  •)  is  a  multi-index  with  non-negative  integer  entries  and  HQfe(x)  is 
the  Hermite  polynomial  of  order  (see  Section  2  for  more  detail).  The  WCE  solution 
of  the  equation  is  given  by  the  series  u  (x)  =  J2aej  (*)  £a,  where  ua  =  E  [«£«]  .  One 
can  view  the  Cameron-Martin  expansion  as  a  Fourier  expansion  that  separates  random 
and  deterministic  components  of  the  equation.  We  demonstrate  that  the  deterministic 
components  ua  (x)  are  uniquely  defined  by  a  lower  triangular  system  of  linear  deterministic 
elliptic  equations  (see  (3.5)  below).  We  refer  to  this  system  as  the  (uncertainty)  propagator. 
Since  the  basis  elements  are  given  by  explicit  formulas,  the  propagator  is  the  key  element 
of  the  solution. 


From  the  numerical  stand  point,  it  is  important  that  the  lower  triangular  (in  fact,  bi¬ 
diagonal)  structure  of  the  propagator  ensures  linear  computational  complexity.  While  it 
is  possible  to  define  an  approximate  WCE  solution  for  the  standard  multiplicative  model 
(1.1),  the  corresponding  propagator  turns  out  a  full  system. 

From  the  statistical  point  of  view,  it  is  important  that  a  solution  of  the  Wick  product  model 

(1.4)  is  an  unbiased  random  perturbation  of  a  solution  of  the  deterministic  equation  (1.5) 

—V  (E  [a  (x,u)]  Vv(x))  =  f(x),  x  e  O  C  Md,  u\qo  =  g  (x)  (1.5) 

in  that  the  mean  (expectation)  of  the  solution  of  equation  (1.4)  is  a  solution  of  equation 

(1.5) . 

In  this  paper  we  investigate  elliptic  equations  of  arbitrary  deterministic  and  stochastic  order. 
As  usual,  by  deterministic  order  we  understand  the  order  of  the  highest  differential  operator; 
in  particular,  equation  (1.1)  has  deterministic  order  two.  The  “stochastic  order”  of  the 
equation  is,  by  definition,  the  highest  order  of  the  Hermite  polynomial  in  the  Wiener  chaos 
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expansion  of  the  random  coefficients.  For  equation  (1.1)  with  the  coefficient  a(x,uj)  defined 
by  (1.3)  the  stochastic  order  is  1;  for  the  same  equation  with  the  log-normal  coefficient  (1.2) 
this  order  is  infinite.  Note  that  any  a(x,uj )  which  is  nonlinear  function  of  £  :=  {£*.  (o;)}fc>1 
has  stochastic  order  higher  than  1. 

The  stochastic  order  is  an  important  and  also  intricate  characteristic  of  the  equations  in 
question.  For  example,  analytically,  problem  (1.1)  with  Gaussian  a(x,u)  (stochastic  order 
1)  is  much  more  difficult  than  the  lognormal  a(x,u>)  (stochastic  order  oo).  One  reason  for 
this  is  that  in  the  Gaussian  case  problem  (1.1)  is  ill  posed  (from  the  classical  point  of  view). 
However,  as  we  will  see  later,  this  ill  posedness  is  irrelevant  in  our  setting.  In  particular,  it 
is  shown  in  Section  3  that  bilinear  stochastic  elliptic  PDEs  are  uniquely  solvable  under  a 
few  simple  conditions.  The  key  is  the  following: 

Assumption  A.  The  expectation  of  the  highest  order  (differential)  operator  is  a  nondegen- 
rate  elliptic  operator. 

For  example,  in  the  setting  (1.1)  this  means  E  \a(x,La)\  >  c  >  0  . 

In  Section  4  it  is  shown  that  under  assumptions  that  are  very  close  to  those  that  guaran¬ 
tee  the  existence  and  uniqueness  of  the  solution  one  can  construct  effective  finite-element 
approximation  algorithms.  More  specifically,  using  the  lower  triangular  structure  and  the 
linearity  of  the  uncertainty  propagator  and  incorporating  the  estimates  of  the  operator 
norms,  we  provide  an  a  priori  error  estimate  for  the  convergence  of  spectral/Iip  finite  ele¬ 
ment  method.  For  equations  of  stochastic  order  one,  this  analysis  has  been  carried  out  in 
[22],  The  same  strategy  can  be  replicated  to  obtain  new  a  priori  error  estimates  for  other 
numerical  methods  for  solving  the  propagator. 

As  mentioned  above  the  WCE  solution  is  given  by  the  series 

u(x)  =  ^2  U<x{x)Ha(€)/Vcd.. 

a&J 

However,  very  simple  examples  (see  Section  3)  demonstrate  that  the  variance  of  WCE 
solutions  is  typically  infinite  or,  more  precisely,  that  E  ||«||2  =  Yla&j  llnc*||2  =  °°-  Fortu¬ 
nately,  the  blow-up  of  the  Wick  product  model  is  controllable,  in  that  the  WCE  solution 
can  be  effectively  rescaled  by  simple  weights  ra  such  that  sf2,OLr2au2a  <  oo  (see  Theorem 
3.3).  We  remark  that,  in  the  “standard  product”  setting,  the  blow-up  of  the  WCE  solution 
is  typically  much  more  severe  than  in  the  Wick  product  setting  (see  Section  5) 
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There  are  alternative  ways  to  address  bi-linear  elliptic  SPDEs.  The  most  developed  al¬ 


ternative  approach  is  based  on  Hida’s  white  noise  analysis  [7].  The  white  noise  approach 


exploits  the  built-in  set  of  stochastic  spaces,  such  as  Hida  or  Kondratiev  spaces  [11,  12], 
or  even  larger  exponential  spaces  [16].  The  traditional  approach  [17,  20,  21,  etc.]  has  to 


select  a  stochastic  space  and  then  to  study  the  largest  possible  class  of  equations  admitting 


a  solution  in  that  space.  The  difference  of  our  approach  is  that  we  select  the  stochastic 
space  that  is  in  some  sense  optimal  for  the  particular  equation  under  consideration. 

The  paper  consists  of  four  sections.  Section  2  describes  the  classes  of  solution  spaces  for 
the  equations  of  interest.  The  theorem  about  existence  and  uniqueness  of  the  solution  is  in 
Section  3.  Numerical  analysis  of  the  equation  using  finite-element  approximation,  including 
the  proof  of  convergence,  is  in  Section  4.  Finally,  Section  5  presents  results  of  numerical 
experiments. 


2.  Weighted  chaos  spaces 


Let  J  be  the  set  of  multi-indices  a  =  (aq,  a  2, . . .)  such  that  a &  E  {0, 1,2,.. .}.  For  a,  (3  E  77, 
we  define 


k>  1 


k>  1 


By  definition,  a.  >  0  if  |ck|  >0  and  (3  <  a  if 


Pk  <  ctk  for  all  k  >  1. 


If  (3  <  cc,  then  we  define 

a  -  (3  =  (aq  -  Pi,  a2  -  /32, . . .). 
We  use  the  following  notations  for  the  special  multi-indices: 


(1)  (0)  is  the  multi-index  with  all  zero  entries:  (0)fc  =  0  for  all  k; 

(2)  e(i )  is  the  multi-index  of  length  1  and  with  the  single  non-zero  entry  at  position  i: 


1,  if  k  =  i : 

e(i)k  =  < 

0,  if  k  ^  i. 

We  also  use  convention  e(0)  =  (0). 


Given  a  sequence  q  =  (%,  k  >  1)  of  positive  numbers  and  a.  E  J ,  we  define 
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Next,  we  introduce  the  following  objects: 

(1)  F  =  (fi,  J7,  P),  a  complete  probability  space  that  is  big  enough  to  support  countably 
many  independent  Gaussian  random  variables; 

(2)  £  =  k  >  1},  i.i.d.  standard  normal  random  variables  on  F; 

(3)  V.  a  real  separable  Hilbert  space; 

(4)  Hn  =  Hn(f),  n  =  0, 1,  2, ,  i  £  R,  one-dimensional  Hermite  polynomial  of  order  n: 


(2.2) 


In  particular,  Ho(£)  =  1,  Hi(f)  =  t,  H2(£)  =  t2  —  1,  H3(f)  =  ti  —  3 1, . . . 

(5)  Stochastic  Hermite  polynomials 


h«(o= 


k>  1 

For  example,  if  a.  =  (0, 2, 0, 1,3,  0,  0, . . .),  then 

Hc(o  =  (e22-i)e4(e53-3e5). 

Recall  the  following  result  (Cameron  and  Martin  [3]). 
Theorem  2.1.  The  collection 


is  an  orthonormal  basis  in  the  space  of  square-integrable  random  variables  that  are  measur¬ 


able  with  respect  to  the  a-algebra  generated  by  the  collection 


In  what  follows,  assume  that  the  cr-algebra  J-  in  F  is  generated  by  £,  and  denote  by  L2(¥]  V ) 
the  collection  of  square-integrable  F-valued  random  elements.  By  Theorem  2.1,  every  v  E 
L2(F;  V )  has  a  unique  representation 


(2.3) 


ot£j 


where 


E  [t>Ha(fl] 

rv! 


(2.4) 


and 


OL&J 


(2.5) 
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Then,  given  a  collection  1Z  =  {ra,  a  G  J"\  of  positive  real  numbers,  we  define  the  space 
1ZL,2(¥'1  V)  as  the  closure  of  L2(¥;  V)  in  the  norm 

IM&L2(F;Vj  =  'Y2  raQ:-  llw“llv-  (2-6) 

a£j 

The  space  1ZL2(¥]  V )  is  called  a  weighted  chaos  space.  We  use  the  notation 

KL2(¥)  =  ^L2(F;M). 

Proposition  2.2.  Let  {ua;  ot  G  J}  be  a  collection  of  elements  from  V.  Then  there  exists 
a  collection  1Z  such  that  YlaeJ  £  TZL2(¥\  V). 

Proof.  We  use  the  following  result:  if  p  >  1,  then 

(2N)pc*  <  °° 

aej  v  ' 

where  N  is  the  sequence  of  positive  integers  and 

(2N)P“  =  JJ(2  k)pak 

k>  1 

(see  [8,  Proposition  7.1]).  Then  it  is  enough  to  take 

1 

r“  “  a!  (2N)2«(1  +  pUjp’ 

□ 

Definition  2.3.  A  generalized  random  element  is  a  collection  {ra,  a  E  J}  of  elements 
from  V. 

3.  Stochastic  Elliptic  Equations:  Existence  and  Uniqueness  of  Solutions 

Let  X  be  a  subset  of  J  such  that  (0)  G  X  and  let  V'  be  a  real  separable  Hilbert  space. 
Consider  a  collection  {Aa,  a:  G  X}  of  bounded  linear  operators  from  V  to  V' ,  such  that  the 
operator  A(0)  has  a  bounded  inverse.  The  objective  of  this  section  is  to  study  the  existence 
and  uniqueness  of  solution  of  the  following  stochastic  equation 

^(A^)oH^)  =  /  (3.1) 

pel 

where  /  G  V'  is  non-random  and  o  denotes  the  Wick  product  [8] .  In  particular, 


Ha«)oH^)  =  Ha+^({). 


(3.2) 
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Definition  3.1.  Equation  (3.1)  is  said  to  have  finite  stochastic  order  N  if 

max  I  a  I  =  N. 

a£l 

The  equation  (3.1)  is  said  to  have  infinite  stochastic  order  if 

sup  |a|  =  oo. 
c*ex 

A  special  case  of  (3.1)  is 

Au  +  y^Mkiiotk  =  /;  (3.3) 

fc>  l 

this  equations  is  studied  in  [14].  In  (3.3),  AQ  =  0  for  |ck|  >  1,  and  the  equation  has 
stochastic  order  1: 

A(o)«  +  E  (Ae(fc)u)oH6(fc)  (£)  =  /.  (3.4) 

fc>  l 

In  general,  (3.1)  can  be  written  as 

A(0 )U  =  f  -  (APU)  °  HM£) 

/3eX,|/3|>0 

and  interpreted  as  a  bi-linear  stochastic  perturbation  of  the  deterministic  equation  A(0)U  = 
/;  this  is  the  reason  for  assuming  that  (0)  £  X.  The  use  of  the  Wick  product  ensures  that 
the  perturbation  has  zero  average. 

In  view  of  Proposition  2.2  and  property  (3.2)  of  the  Hermite  polynomials,  the  following 
definition  is  natural. 


Definition  3.2.  A  solution  of  equation  (3.1)  is  a  generalized  random  element  such  that 

A(o)m(o)  =  /,  | ck|  =  0, 


A(0 )Ua  =  -  ^2  A/3ua—(3,  M  >  o. 


(3.5) 


/3ex 
0  <(3<a 


Indeed,  writing  u  =  uq;HCi(^),  substituting  formally  into  (3.1),  and  using  (3.2),  we 

find 

^  ^  y  ApUa Hot_[_jg(^)  /, 

ol&j  /3ex 


or 


/ 


E 

aej 


\ 


J2  apu 


ol—(3 


/3ex 
\/3<« 


H«(£)  =  f- 


Equating  the  coefficients  of  every  Ha(^),  we  conclude  that  (3.5)  must  hold. 
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The  following  theorem  establishes  the  existence  and  uniqueness  of  solution  of  equation  (3.1). 
Theorem  3.3.  Assume  that 

(1)  The  operator  A(0)  has  a  bounded  inverse  from  V'  to  V. 

(2)  Each  of  the  operators  A^A^,  (3  €  1,  \/3\  >0,  is  bounded  on  V,  with  the  operator 
norm  Cp. 

(3)  There  exists  a  sequence  b  =  {b\,  62 1  •  •  •)  of  positive  numbers,  such  that,  for  all  (3  £  T, 

Cp  <  bP.  (3.6) 

(4)  The  function  f  is  deterministic  and  belongs  to  V' . 

Then 

(1)  equation  (3.1)  has  a  unique  solution  u; 

(2)  the  solution  is  an  element  of  the  weighted  chaos  space  1ZL2( IF;  V)  with 

ca 

ra  =  — T 

a! 

for  some  sequence  c=  (ci,  C2, . . .)  of  positive  numbers. 


Proof.  To  begin,  let  us  first  understand  the  structure  of  the  system  of  equations  (3.5). 
Define 

Bp  =  -A^Ap,  (3  6l 

Then  (3.5)  becomes 

«(o)=A(01)/,  |a|  =  0, 

ua=  ^  B pua-p,  |a|  >  0.  (3-7) 

/3ex 
0  <(3<a 

Thus,  ua  with  |ct|  =  n  >  0  is  determined  by  ua  with  |a|  <  n,  which  implies  existence  and 
uniqueness  of  solution  of  (3.5),  and  hence  of  (3.1). 


In  view  of  (2.7),  to  prove  that  the  solution  is  an  element  of  a  weighted  chaos  space  with 
weights  ra  of  the  form  ra  =  ca /a\  for  some  sequence  c  =  (ci,C2, . . .)  of  positive  numbers, 
it  is  enough  to  show  that  there  exists  a  sequence  of  positive  numbers  q  =  (qk,  k  >  1)  such 
that 


Il«t»||v  <  9cl||'W(0)l|u;  (3-8) 

then  we  can  take  Ck  =  (2 kqk)~2  Thus,  to  complete  the  proof,  it  remain  to  establish  (3.8). 
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We  start  with  the  particular  case  (3.3).  Even  though  this  equation  has  been  studied  in  [14], 
the  analysis  there  relies  on  a  closed-form  expression  for  ua  and  does  not  easily  extend  to 
higher  orders.  Denote  by  B*.  the  operator  A-1Mfc.  Then,  for  | ck |  >  0,  (3.7)  becomes 

Ua  'y  ^  e(k)  ■  (3-9) 

k>  1 

We  now  show  by  induction  that  (3.8)  holds  with 

qk  =  Ak2bk. 


Since  ue (&)  =  BfcU(0)>  we  have  (3.8)  for  |o;|  =  1.  Now  assume  that  (3.8)  holds  for  all  a  with 
|  ck |  <  n  —  1.  Then  (3.9)  and  the  triangle  inequality  imply,  for  a  with  |ck|  =  n, 

jkallv  <  II-(o)IIyEW°  hq01-^ 

k>  1 


qaWu(o)\\v  j 


1  00  i 

ELak>0  ^  „cL\\nt  ||  1  / 


k>  1 


k=  1 


4  k2 


q  IIu(o)IIy, 


because  ^fc>1  k  2  =  (tt2/6)  <  4.  This  completes  the  proof  of  (3.8)  for  equation  (3.3). 


Let  us  now  consider  the  general  equation  (3.1).  We  will  use  assumption  (3.6)  to  show  that 
(3.8)  holds  for  N  >  1  with 

qk  =  C0k2bk, 


where 

c»  = 2  E  (3-io) 

a£j  V  ' 

see  (2.7).  For  the  curious  reader,  we  note  that  numerical  computations  show  that  Co  ~  3.2, 
and  simple  analytic  considerations  put  the  value  of  Co  in  the  range  [3.0,  3.4], 


As  in  the  case  of  equation  (3.3),  the  proof  of  (3.8)  is  by  induction  on  |o;|.  Let  ko  =  min{|/3|  : 
/ 3  £  1}.  If  |q|  <  ko,  then  ua  =  0  and  so  (3.8)  trivially  holds.  If  a  £  1  and  |a|  =  ko,  then 
(3.5)  implies 

Ua  =  Ba'U(o),  ||na||  <  b“||rt(0) ||v  <  qa\\u(o) || v- 

Now  assume  that  (3.8)  holds  for  all  ua  with  |a|  <  n—  1  and  fix  |cc|  =  n  >  ko-  By  (3.7),  the 
triangle  inequality,  and  the  induction  assumption, 

qa  j 

||««||y<  E  IIB/3“«-^llv  <  lho)II^E  (2mw 

/3SJ  /3ex v  ’ 

0</3  <a 

<  ll«(0)ll E  (2N)2/3  - 
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where  the  last  inequality  follows  from  (3.10).  This  completes  the  proof  of  Theorem  3.3.  □ 

Below,  we  illustrate  Theorem  3.3  on  several  examples.  To  emphasize  the  main  points,  we 
consider  a  simplified  version  of  equation  (3.1),  in  which 

(1)  V  =  V'  =  M,  so  that  every  operator  is  multiplication  by  a  real  number;  then, 
with  no  loss  of  generality,  we  take  A(0)  =  1. 

(2)  /  =  I- 

(3)  the  operators  A^  are  such  that  A^  /  0  only  if  (3  =  (n,  0,0,.. .)  for  some  n  = 

0,1,2, . 

In  other  words,  (3.1)  becomes 

u  =  1  +  ^afenoHfc(£),  (3.11) 

k>  1 

where  £  is  a  standard  Gaussian  random  variable.  The  solution  u  has  the  form 

u  =  y~\nH  n(£) 

n>  o 

for  some  real  numbers  un,  and 

E|tt|2  =  ^  n!  (3-12) 

n>  1 

As  a  first  example,  consider  the  equation 

u  =  1  +  mo  £.  (3.13) 

It  is  easy  to  verify  that  u  =  1  +  X^fc>i  Hfc(£)  is  a  solution,  that  is,  Uk  =  1  for  all  k.  By 
(3.12),  u  i  L2(¥]R). 

The  coefficients  of  the  solution  of  (3.13)  are  bounded.  With  a  simple  modification  of  (3.13), 
we  can  get  coefficients  that  grow  exponentially  fast.  Indeed,  consider 

u  =  1  +  uo  £  +  u  o  (£2  —  1);  (3-14) 

recall  that  £2  —  1  =  H2(£).  Then  (3.5)  implies 

V>§  —  1,  U\  —  1,  Un  —  Wn- 1  +  Un— 2,  2. 

Thus,  un  is  nth  Fibonacci  number,  and,  as  n  — ►  oo,  un  behaves  asymptotically  as  (1  + 
\/5)n/2n.  In  other  words,  we  have  un  >  rn  for  some  r  >  1.  Again,  (3.12)  implies  IEu2  =  oo 
and  so  u  ^  L2  (F;  M) . 
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Examples  (3.13)  and  (3.14)  show  that,  if  equation  (3.1)  has  finite  stochastic  order,  then,  in 
general,  the  solution  does  not  belong  to  L2(F;  V).  In  fact,  we  are  not  aware  of  any  elliptic 
equation  with  multiplicative  finite-order  noise  and  with  a  square-integrable  solution. 

Let  us  consider  now  the  equation 


u  =  i  — 


u  o  Hn(£) 


n>  1 


n\ 


(3.15) 


This  equation  has  infinite  stochastic  order.  Direct  computations  show  that 


u  = 


i  +  Etx:H»K)  =  e_e_(1/2>' 


n>  1 


n 


one  can  also  derive  this  by  re-writing  (3.15)  as 

e£-(V2)  o  u  =  i 


and  noticing  that  d/2)  oe  ^  d/2)  =  1.  In  particular,  E u2  =  e  <  oo,  so  u  £  L2(F;M). 


4.  Numerical  Analysis  of  a  Stochastic  Dirichlet  Problem 


Let  O  E  be  a  bounded,  open,  connected  set  with  a  Lipschitz  continuous  boundary  dO. 
Denote  by  (•,  •)  the  inner  product  in  L2(0).  Consider  the  following  problem: 


A  u  +  ^  A/3moH^)  =  /,  xeO, 

/3ei 

I/3|>1 

u  =  0,  x  €  dO , 


(4.1) 


where  £  is  a  finite-  or  infinite-dimensional  standard  Gaussian  random  variable.  For  simplic¬ 
ity,  we  assume  that  the  force  term  f(x)  is  a  deterministic  function.  The  operators  A  =  A(0) 
and  Ap  are  defined  by 


A u(x)  =  -  ^2  Di{aij(x)Dju(x)), 
i,j= 1 
d 

A p{x)  =  wp(x)  Di(aij}p(x)Dju(x)), 
i,j= 1 


(4.2) 

(4.3) 


where  the  functions  aij(x),  aijtp(x),  and  wp(x)  satisfy  the  following  assumptions: 

Assumption  4.1. 


(a)  The  functions  aij(x )  are  measurable  and  bounded  in  the  closure  O  of  O; 
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(b)  There  exist  positive  numbers  A\,  A2  such  that  Ai|y|2  <  aij(x)yiyj  <  A2\y\2  for  all 
x  £  O  and  y  £  Md; 

(c)  The  functions  a^pix)  are  bounded  and  measurable: 

Wij,p(x)  I  <  cp  for  all  x  £  6,  i,j  =  1, . . . ,  d]  (4.4) 

(d)  T/ie  functions  wp(x)  are  bounded  and  Lipschitz  continuous: 

\wp{x)\<Cp  and  \wp(x)  -  wp(y)\ <  Cp\x  -  y\,  x,yeO.  (4.5) 

(e)  There  exists  a  sequence  b  =  (b\, 62  >  •  •  •)  of  positive  numbers  such  that,  for  all  (3  £  Z, 

{Cp  +  Cp)Cp  <  bP. 


Taking  as  V  and  V'  the  Sobolev  spaces  Hq(0)  and  i4_1((9),  respectively,  we  note  that 
Theorem  3.3  applies  to  equation  (4.1)  because  the  operators  A~1Ajg  are  bounded  on  V. 
Indeed,  note  that 

\\u\\v  =  IMll2(e>)  +  l|V«||i2(0) 

and  let  CV  be  the  Poincare  constant  for  the  domain  O: 

IHI L2(0)  ^  CWI|Vt*||ia(0).  (4.6) 


Assumption  4.1(b)  implies 


Ai 

(A u,u)  >  \\u\\v,  ueV- 

1  i  Otj- 


here  (•,  •)  denotes  the  duality  paring  between  V'  and  V  relative  to  the  inner  product  in 
L2(0),  Then,  by  direct  computation, 


II A  1Apu\\v  < 


1  +  C-„ 


(C%  +  C%)C£. 


Example  4.2.  Let  G(x)  is  a  Gaussian  field  on  O  with  representation 

OO 

G(x)  =  ^2hk(x)£k, 
k= 1 

where  {hk,  k  >  1}  are  smooth  bounded  functions  such  that 

Y  (  sup  \hk(x)\2  +  sup  \Vhk{x)\2  )  <  00 
Y\  xeO  J 


(4.7) 


(4.8) 


and  £  =  {£fc,  k  >  1}  is  a  collection  of  independent  standard  Gaussian  random  variables. 
Note  that 

1 


G(x)  :=  EeG(a:)  =  exp  -  Y  h2k{x) 


k>  1 
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Examples  of  (4.7)  are  Karhuan-Loeve  expansion  [13],  spectral  expansion  [19],  and  Fourier- 
Wavelet  expansion  [6]. 

Consider  the  equation 

V  •  (eG(a,)  o  Vu)  =  /.  (4.9) 


If 

then 


M*)  =  “rll^C®)’  a<E^’ 

fc 


eG^  =  G(®)(l  +  M*)H/3(^))- 

I/3|>1 

and  (4.9)  becomes  a  particular  case  of  (4.1): 


V  •  (G(x)Vu)  +  V  •  (G(x)f}P(x)Vu)  oH^(|)  =  /; 
)/3|>l 

to  satisfy  Assumption  4.1(e),  take  b k  =  sup^gg,  \hk(x)\  +  sup^g^  \\7hk{x)\. 
Replacing  in  (4.9)  the  usual  exponential  eG ^  with  the  Wick  exponential 

eoG((i))  _  eG(x)-±J2k>ihl(x) 

also  leads  to  an  equation  of  the  form  (4.1). 


(4.10) 


The  objective  of  this  section  is  to  study  numerical  solution  of  (4.1).  Specifically,  we  construct 
an  approximation  of  (4.1)  using  finite  elements  in  the  physical  space  O  and  a  truncated 
chaos  expansion  in  the  probability  space.  For  the  sake  of  concreteness  we  assume  that  the 
physical  space  is  two-dimensional:  d  =  2. 


Let  ^  be  a  family  of  triangulations  of  O  with  straight  edges  and  h  the  maximum  size  of 
the  elements  in  We  assume  that  the  family  is  regular,  in  other  words,  the  minimal 
angle  of  all  the  triangles  is  bounded  from  below  by  a  positive  constant.  We  define  the  finite 
element  space  as 


VhK  = 


[v  |  V  o  Ff)1  e  &p(R)},  vh,p  =  [v  e  H\D)  |  v\k  g  VhKp,  K  e  .%), 


where  Fk  is  the  mapping  function  for  the  element  K  which  maps  the  reference  element 
R  (for  example,  an  equilateral  triangle  or  an  isosceles  right  triangle)  to  the  element  K 
and  £?p(R)  denotes  the  set  of  polynomials  of  degree  at  most  p  on  R.  We  assume  that 
Vh\dO  =  0,  Vh  £  Vh,p.  Thus,  Vh,p  is  an  approximation  of  Hq(0)  by  piece-wise  polynomial 
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functions.  There  exist  many  choices  of  basis  functions  on  the  reference  elements,  such  as 
h- type  finite  elements  [5],  spectral//ip  elements  [10,  18],  etc. 

To  describe  truncated  chaos  expansion,  fix  positive  integers  M,  n  and  define  the  finite  di¬ 
mensional  set  JM,n  C  J  as 


Jm,u  ■=  {«  H  <  M,  a.k  =  0,  k  >  n|. 
The  truncated  Wiener  chaos  space  VM,n  is  then 

VM’n:={  J2  faPM  |  /«GR}, 

Similarly,  let 

-^M,n  —  2^  n  • 


(4.11) 


(4.12) 


It  follows  from  (3.5)  that  ua  is  fully  determined  by  coefficients  up  and  operators  Ap  with 
/3  <  a.  Thus  we  only  need  to  keep  the  operators  A p  with  /3  G  Zm.u  if  we  only  consider  the 
coefficients  ua  with  a  G  Jm,u- 

The  sFEM  (stochastic  finite  element  method)  can  be  formulated  as  follows:  Find  € 
VhtP  0  VM’n  such  that,  for  all  v  G  Vh,p  0  VM,n, 


E(A(o .)«££» +E  £  (Apu%?oHp{Z),v)  =E(f,v).  (4.13) 

|/3|>0 

Remark  4.3.  The  underlying  idea  of  sFEM  is  to  consider  a  finite  element  approximation, 
followed  by  a  truncation  of  the  chaos  expansion.  Alternatively,  it  is  possible  to  truncate  the 
chaos  expansion  first  and  then  apply  the  finite  element  method  to  the  truncated  version  of 
the  uncertainty  propagator  (3.5): 


A(o)«(o)  =  /■  |«|  =  o, 

A(o)^c*  +  ^  Apua-.p  =  0,  |o:|  >  0.  (4.14) 

0  <(3<a 

This  approach  will  not  be  discussed  in  this  paper. 


The  following  is  the  main  result  about  the  convergence  of  the  stochastic  finite  element 
method. 
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Theorem  4.4.  Assume  that  conditions  given  in  Assumption  4-1  and  Theorem  3.3  are 
satisfied  and  in  addition,  u  E  1ZL2(¥]  Hm+1(0))  for  some  integer  m  >  1.  Then 

II  M,n  1 1  ^  7  u —  1  —mil  || 

\\u  —  Uh,p  \\kL2(F-,H£(0))  —  P  \\u\\'R.L2(¥-)Hm+1(0)) 


+  C2 1| /||  77-1(0)  (  p?  + 


(4.15) 


where  p  =  min(m+  1  ,p  +  1),  the  constants  C\  and  C2  are  independent  of  h  and  p,  and  C2 
is  also  independent  of  M,  n. 


Proof.  Define  V  =  Hq(0)  and  let  the  Wiener-chaos  expansion  of  u  and  be 

u=  ^2  UoMofiOi  uh,pn  =  ^2  “a Ha(£), 


Ot&J 

M,n 

respectively.  Then 

M.n  \  '  t 

u~uh,p  =  2s  (Ua~ 

OteJ\Ju,n 

and 

II  M,nii2 

1  \u~uh,P  \\tzl2{¥-v)  ~  2s 

II  ^  1 1 2  1 1 TT  1 1 2 

\\uol  —  uol\\v  lln«II^L2(F) 

+  E 

uc* \\v  Pa  7?X2(F)  =  h  +  I2 

ateJ\JM,n 


The  structure  of  equation  (3.5)  ensures  that  the  coefficients  ua  with  a  E  J\Jm,u  and 
operators  A^g  with  (3  E  T\Zj/[,n  do  not  affect  ua  with  ol  £  Jm,u-  In  other  words,  I\  is 
the  error  from  the  finite  element  discretization  and  I2  is  from  the  truncation  of  the  chaos 
expansion. 

We  first  look  at  the  term  I\  by  considering  equation  (4.14).  Since  the  equation  for  ix(0) 
is  a  usual  deterministic  elliptic  equation,  we  here  only  consider  the  equations  for  ua  with 
|a|  >  0.  By  (4.13), 

A(o  )Ua+  ^2  A^cx-f3  =  0.  (4.18) 

f3£^M,n 
0</3  <ct 

Furthermore,  although  we  truncate  the  number  of  random  variables  and  the  multi-index 
set  Z,  the  structure  of  the  uncertainty  propagator  (3.5)  ensures  that  the  exact  solution  also 
satisfied 

A(o)^ct  T  ^  '  A f^ua—f3  0. 

0</3  <a 


(4.19) 
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Denote  by  a(-,  ■)  the  bilinear  form  generated  by  the  elliptic  operator  A(0): 


a(u,  v )  =  (A u,  v). 


Then,  for  every  G  V/,, 


Q>(ua  VfnUa  Vfi ) 

=  a(ua- vh,ua~  vh)  +  a(ua- ua,ua- vh) 

=  a(ua  -  vh,  ua  -  vh)  -  ^2  (ap(u<x-(3  -  ~  Di>  •  (4.20) 


/3exM,n 

0<f3<a 


Assumption  4.1(b)  implies 


A 


i  +  a 


tHMIu  <  a(u,u),  a(u,v)  <  A2||rt||vHM|y,  u  G  V, 


where  Cn  is  the  Poincare  constant  for  the  domain  O;  see  (4.6).  Then,  with  C\  =  A\/(\+Cn) 
and  thinking  Wh  =  ua  —  Vh, 


Ci\\ua  -  vh\\v  <  A2\\ua  -  vh\\v  +  sup 

wh^vh,P 


^  '  (Ajg(?iQ:_^g  p) ,  Wh) 


0</3<a 


V 


Combining  the  above  inequality  with  the  triangle  inequality 


|  Ua  ~  Ua\\v  <  1 1 He*  —  Vh  ||  V  +  H^ct  ~  Vh\\v  i 


we  find 


|  Vot  Uct  1 1 U  A 


1  +  7W  )  inf  \\ua-VhWv 

c  i  J  vhevh }P 


1 

+  —  sup 
Ul  wheVKp 


^  '  (A f3(ua—/3  Hct— ft) ,  Wh) 


0  </3<a 


V 


(4.21) 


By  Assumptions  4.1(c,d) 
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where  Cp  =  (7g((7g  +  Cp) .  Therefore, 

\\ua  Wa||v  E  C  inf  || U(x  n/i 1 1  w  T  ^  '  Cp\\ua—p  Wet— /3||v 

Vh  £Vh,v 

0</3<« 


(4.23) 


for  some  constants  C  and  Cp.  The  first  term  on  the  right-hand  side  is  the  approximation 
error  of  finite  element  approximation  and  the  second  term  is  the  error  propagated  from  the 
approximation  of  lower-order  coefficients.  Using  equation  (4.23)  recursively,  we  obtain 


\ua  ~  ua\\v  <  Y  Catp  inf  \\up-vh\\v, 

z - -  Vh&Vh  v 

(3<cx  ,p 


(4.24) 


with  some  constants  Cap.  According  to  the  theory  of  the  finite  element  approximation 
[1,  5],  there  exists  a  positive  number  Cf  such  that,  for  every  function  u  £  Hq  (0)r\Hm+1(0), 


inf  \\ua-Vh\\ v<Cfh^  Lp  m\\ua\\Hm+i(D), 
vhevh  v  1 


(4.25) 


where  p  =  min (m  +  1  ,p  +  1).  Then 


\ua-UocWv  <  C2h2fl  2p 


—  clnr—  2m 


< 


^  ^  C QL,f3  1 1  ^OL 1 1  fjrn+l  (£)) 
(3<CX 

Cjh^-2p-2m  £  dip  Y,  \M2Hm+HO) 

/3<a.  /3<c* 


—  k 2/i  2p  2mCa  llU/3||/fm+l(<9), 


(3<cx 


and  the  term  I\  can  be  bounded  as 


h  <  h2fl-2p-2m 


< 


E  1 1  UP  1 1  Irn  +  l  (£))  1 1  H„  1 1  lL2  (F) 

\J M,n  ft— CX 

h2»-2p-2mC(M,n)  Y  ^\\ua\\jjrn+l{D)ma\\lL2{¥) 

CxdzU M,n 


^  h  ^  p  C(Mj  n)  max  Col  ^  ^  ll^«lli/m+1(L))ll^'Q:ll7^i/2(F) 

,n  _ 

Ol(E:u  M,n 

<  h2^  2p  2mCi(M,  naMI^L2(F;tfm+i(D))- 


(4.26) 


To  estimate  the  truncation  error  I2  of  the  Wiener  chaos  expansion,  note  that  ||Hce||^£2^ 
ra  cd,  and,  by  (3.8), 

 (2N)“2o:g-Q 


«a||y  <  <?°1|w(o)||y  <  C'Aq^WfWH-^iD)- 


Also 
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As  a  result, 

h=  Y.  IKIlv'IIHJIk.m  <  CxII/IIh-i(d)  E  (2N)-2“<c2l^  +  i), 

«Si7\v7M,n 

where  the  last  inequality  is  taken  from  [4,  proof  of  Theorem  2],  This,  together  with  (4.17) 
and  (4.26),  completes  the  proof  of  Theorem  4.4.  □ 


5.  Computations 

In  this  section  we  present  some  results  of  numerical  experiments  and  compare  our  model 
with  another  one  studied  in  the  literature  on  uncertainty  quantification.  We  consider  the 
following  one-dimensional  problem 

)  dx\  dx  >  ’  (5.1) 

[  u(0)  =  0,  u(l)  =  0, 

where 


•  *  denotes  either  the  Wick  product  o  or  the  ordinary  product  •  ; 

•  K  =  ec^-2c“,  £  is  a  standard  Gaussian  random  variable,  c  >  0  is  a  real  number. 

The  random  variable  K  can  be  regarded  as  a  simplified  version  of  the  random  field  eoG^ 
from  (4.10).  Note  that 

E  [K]  =  1,  Var  [K]  =  -  1.  (5.2) 


Denote  by  Hj  the  normalized  Hermite  polynomials: 

m) 


Hi(0  = 


Substituting  the  truncated  Wiener  chaos  expansion 

n 

u(x)  =  ^Ui(x)Hj(£) 

i=0 

into  equation  (5.1)  and  implementing  the  finite  element  projection  we  obtain  the  uncertainty 
propagators 

n  ,2 

~  ^  \KRjRj]  =  Sqj,  j  =  0,1, . . .  ,n  (5.3a) 

i=o  a‘ 

n  d? 

-  ^2  -j-^u.i(x)E  [(K  o  H,;)Hj]  =  S0j,  j  =  0, 1, . . . ,  n  (5.3b) 

i=o  ax 

for  the  ordinary  product  and  the  Wick  product,  respectively. 
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Using  the  Taylor  expansion  and  the  properties  of  Hermite  polynomials,  we  obtain 


oo  [n/2] 

E  [A'EjHj]  =  Y,  E  E  A(n,  k)B(i,j,p)5n-2k,i+j-2P,  (5.4a) 

n= 0  fc=0  p<iAj 


__  _  p-c2/ 2  °°  I”/2]  //  n  k  +  iV 

E  [(A  .  H,)Hj]  =-^EE 

72 — 0  k —  0 


(5.4b) 


where 


A(n,k)  =  —  f  2k  ){2k  -  1)!! y/(n  -  2k)\,  B(i,j,p )  = 


Mr ) (i+j~ 2p 

p)  Xp)\  i~P 


1/2 


Formulas  similar  to  (5.4a)  and  (5.4b)  can  also  be  derived  when 


K{x )  =  eoG^x\ 


(5.5) 


with  G(x)  given  by  (4.7).  For  simplicity,  we  present  numerical  results  for  the  space- 
independent  case  G(x)  =  c£,  corresponding  to  K  =  //2^.  The  qualitative  behavior 

of  the  solutions  of  (5.1),  with  either  ordinary  or  Wick  product,  remains  the  same  for  the 
more  general  noise  (5.5),  as  long  as  conditions  (4.8)  hold. 

Let  us  compare  the  matrix  structures  of  E  [iFHjH/]  and  E  [(A'  o  Hj)Hj] . 

•  Ordinary  product.  For  all  ci, . . . ,  cn  £  M,  we  have 


CiCjE  [ATHjHj] 

i,j= 1 


E 


E 


CiCjKHiHj 

ij= 1 

/  72  N 

V  i=l  2 


>  0. 


It  can  be  shown  that  when  A"  is  lognormal,  c,  must  be  zero  if  the  equality  is  satisfied. 
In  other  words,  the  matrix  E  [ATH,Hj]  is  symmetric  and  positive  definite.  Thus  all 
the  eigenvalues  of  the  matrix  (E  [ATHjHj]  ,  i ,  j  =  1 ,  are  positive,  and  the 

system  (5.3a)  is  well-posed.  However,  all  unknowns  in  equation  (5.3a)  are  coupled 
together,  which  is  not  a  desirable  property  for  numerical  computation. 
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Wick  product.  As  n  must  be  even  for  i  =  j  to  get  nonzero  terms  in  equation 
(5.4b),  it  follows  that  IE  [(/\  o  Hj)Hj]  =  1  for  all  i  =  0, 1, . . ..  Indeed, 

E  [(K  oH,)Hj  =  ^ 52m.2k+i  i 

LV  *}  ^  (2m  -  2k)\(2k)\\  +  ’ 

e-i/2  “  \fi\ 

s[i\  ^  (2m)!! 

v  m=0  x  ' 


=  e-1/2£ 

m= 0 


i 


2mm\ 


=  e-V2ei/2  =  L 


Furthermore,  we  note  that  E  [(A'  o  Hj)Hj]  is  a  lower-triangular  matrix.  In  other 
words,  the  coefficient  Ui  depends  only  on  Uj  with  j  <  i. 


In  the  first  numerical  experiment,  we  take  c  =  1  and  compare  the  chaos  coefficients  of 
solutions  corresponding  to  the  ordinary  and  Wick  products.  Figure  1  presents  the  results. 
For  a  better  comparison,  we  use  the  same  scale  on  both  graphs  and  also  plot  the  function 
Uc].  the  solution  of 

- d =  1’0<x<1’’  ^(°)  =  W)  =  (5-6) 

Note  that  Ud  is  the  solution  of  (5.1)  when  K  is  replaced  by  its  mean  value  E  [K]  =  1. 

We  make  the  following  observations: 


(1)  When  the  ordinary  product  is  used,  Ud  7^  rio!  when  the  Wick  product  is  used, 
Ud  =  no-  This  is  a  direct  consequence  of  equations  (5.3a)  and  (5.3b).  In  Figure  2  , 
we  plot  the  function  uq  corresponding  to  the  ordinary  product,  for  different  values 
of  c.  It  appears  that  u$  is  an  increasing  function  of  c  and  is  close  to  Ud  when  c  <  0.1. 

(2)  The  larger  the  number  i,  the  smaller  the  coefficient  tij.  However,  it  appears  that 
the  coefficients  for  the  Wick  product  are  much  smaller  than  the  corresponding  co¬ 
efficients  for  the  ordinary  product.  Below,  we  investigate  this  phenomenon  further 
by  changing  the  value  of  c. 


In  the  second  numerical  experiment,  we  study  the  relation  between  the  variance  of  the 
perturbation  K  and  the  variance  of  the  corresponding  solution.  We  change  the  value  of  c 
from  0.1  to  2.5,  which,  by  (5.2),  is  equivalent  to  changing  the  standard  deviation  of  K  from 
0.1003  to  22.7379.  Figures  3  and  4  present  the  results. 

In  Figure  3,  we  plot  the  variance  of  the  solution  as  a  function  of  x  for  several  values  of 
c,  and  notice  that,  in  each  case,  the  variance  seems  to  be  largest  near  the  middle  of  the 
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Figure  1.  Chaos  expansion  coefficients.  Left:  ordinary  product;  Right:  Wick  product. 


Figure  2.  The  plot  of  E  [u(x)]  for  different  values  of  c,  ordinary  product. 

interval  (0, 1).  In  Figure  4,  we  plot  the  solution  variance  at  point  x  =  0.5  as  a  function  of  c. 
It  appears  that  variances  corresponding  to  ordinary  and  Wick  products  are  close  for  small 
c,  but  diverge  quickly  as  c  increases:  in  the  same  range  of  the  values  of  c,  the  variance  of 
the  solution  increases  from  O(10~'5)  to  O(106)  for  the  ordinary  product,  and  from  O(10~5) 
to  O(IC)1)  for  the  Wick  product. 

In  summary,  although  the  solutions  of  equation  (5.1)  with  Wick  and  ordinary  products  be¬ 
long  to  L2(E;  iLg(0, 1)),  the  two  different  products  lead  to  many  differences  in  the  properties 
of  the  solution.  The  differences  are  especially  pronounced  when  the  noise  I\  is  large.  A 
more  detailed  analysis  shows  that  many  numerical  characteristics  of  the  solutions  differ  by 
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Figure  3.  Variance  of  the  solution  for  different  values  of  c. 
Left:  Ordinary  product;  Right:  Wick  product 


Figure  4.  Variance  of  the  solution  at  x  =  0.5,  as  a  function  of  c. 

the  factor  of  e°2 ,  which  is  related  to  the  equality  ec^-2c2  o  e_c^_2c~  =  1.  At  this  point, 

2 

we  make  the  statement  intentionally  vague,  because  this  scaling  by  the  factor  ec  seems 
to  have  many  important  ramifications  and  must  be  studied  further,  both  analytically  and 
numerically. 
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